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ELASTOHYDRODYNAMIC LUBRICATION THEORY* 


This chapter is devoted to the method used by the authors 
to solve the elastohydrodynamic lubrication problem for elliptical 
contacts (Hamrock and Dowson, 1976a). The problem has been dis- 
cussed earlier in some detail and the relevant effects examined 
separately, so it remains only to seek a solution that reconciles 
the hydrodynamic and elastic equations. The problem is to calculate 
the pressure distribution in the contact and at the same time to 
allow for the effects that this pressure will have on the properties 
of the fluid and on the geometry of the elastic solids. The solution 
will also provide the shape of the lubricant film, particularly the 
minimum clearance between the solids. Thermal effects are, however, 
neglected. The variation of lubricant viscosity with temperature and 
the thermal expansion of the lubricant and the solids are not 
considered. 

7.1 Reynolds Equation 

The coordinate system to be used in defining the elastohy- 
drodynamic problem is shown in Figure 7.1. The appropriate 


*Published as Chapter 7 in Ball Bearing Lubrication by 
Bernard 0. Hamrock and Duncan Dowson, John Wiley & Sons, Inc., 
Sept. 1981. 





Reynolds equation for the coordinate system (x,y) can be written 
from equation (5.43) as 
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(5.43) 


where u = (u^ + Uj^)/2 is the mean surface velocity or the 
entraining velocity in the x direction. Letting 



(7.1) 


(7.2) 


where U = udq/E'Rj^ is a dimensionless speed parameter and 
k = 1.03 (R /R )^’®^ is an ellipticity parameter. 

j X 

Equation (7.2) is a nonhomogeneous partial differential 
equation in two variables. This is generally a difficult equa- 
tion to solve, and the degree of complexity depends on the form 
of the parametric functions and H and on the boundary 

conditons. In elastohydrodynamic lubrication the parameters are 
of the most difficult form. Therefore equation (7.2) is solved 
by approximate numerical methods. Before proceeding, however, 
it is necessary to write expressions for the dimensionless film 
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thickness H, the dimensionless density p, and the dimension- 
less viscosity "n. 

The influence of pressure on the density of fluids has been 
discussed in Section 5.5. From this section the dimensionless 
density can be written as 


A^PE’ 

p = 1 + r 

1 + B 2 PE’ 


(5.57) 


where and B 2 are constants that depend only on the 
fluid considered. 

The viscosity of fluids has been discussed in Section 5.4, 
where attention was drawn to the work of Roelands (1966). For 
isothermal conditions the dimensionless viscosity of the fluid 
can be written from Section 5.4 as 


n = 



l-(l+p/c) 



(5.52) 


where 

n„= 6.31x10"^ N s/m^ (9.15x10“^ Ibf s/in^) 
c = 1.96x10^ N/m^ (28 440 Ibf/in^) 
riQ = viscosity of lubricant at atmospheric pressure 
= viscosity-pressure index, a dimensionless constant 
These dimensions should also be used to define the constants in 

equation (5.52). 
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7.2 Film Shape 


The film thickness can be written simply as 

h(x,y) = hp + S(x,y) + 6(x,y) (7.3) 


where 

hQ = constant 

S(x,y) = separation due to geometry of undeformed ellipsoidal- 
sol ids 

6(^»y) = elastic deformation 

The separation due to the geometry of the two undeformed 
ellipsoidal solids shown in Figure 2.18 can be described by an 
equivalent ellipsoidal solid near a plane. The geometrical re- 
quirement is that the separation of the two ellipsoidal solids 
in the initial and equivalent situations should be the same at 
equal values of x and y* Tnerefore from Figure 7.1 the sep- 
aration due to the undeformed geometry of the two ellipsoids can 
be written as 


S(x,y) 


(x - mb)^ ^ (v - la)^ 
2Rx 2Ry 


(7.4) 


where 

m = constant used to determine length of inlet region 
a = constant used to determine width of side-leakage region 
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For illustration the mesh described in Figure 7.1 is used, 
where m = 4 and *• = 2. However, the equations developed are 
written in general form. 

Figure 7.2 illustrates the film thickness and its compo- 
nents for an ellipsoidal solid near a plane. Substituting 
equation (7.4) into equation (7.3) while using equation (7.1) to 
make this equation dimensionless gives 



where Hq is a constant that is initially estimated. 

The way in which the elastic deformation of an equivalent 
ellipsoidal solid near a plane can be calculated has been out- 
lined in Section 5.7. Therefore from Figure 7.1 and the results 
of Section 5.7 the elastic deformation can be written as 

_2AXc (m+-n)Xd 

- i 

j— 1» 2 , . • • i=l , 2 , . • • 

where 

number of divisions in semimajor axis 
TT = constant used to determine length of outlet region 
F = number of divisions in semiminor axis 
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i* = |a - i| + 1 

j* = le - j| + 1 



(7.7) 


Equation (7.8) and Figure 7.1 clarify the meaning of equa- 
tion (7.6). The elastic deformation at the center of the rec- 
tangular area 6g ^ (Figure 7.1) caused by the pressure on the 
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various rectangular areas in and around the contact ellipse can 
be written as 



(7.8) 

7.3 Phi (0*) Solution 

Having defined the density, viscosity, and film thickness, 
we can return to the problem of solving for the pressure in the 
Reynolds equation. It is well known (e.g., Whomes, 1966) that 
the dimensionless pressure P of the Reynolds equation plotted 
as a function of x exhibits a very localized pressure field 

p p 

with high values of aP/aX and a P/aX . Such a condition 
with high gradients is not welcomed when performing numerical 
analysis by relaxation methods. Therefore, to produce a much 
gentler curve, a parameter a>* is introduced where 

(7.9) 
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The pressure P is small at large values of film thickness H 
and large when the film thickness is small. This substitution 
also has the advantage that it eliminates all terms containing 
derivatives of the products of H and P or H and «>*. 
Therefore from equation (7.9) and expansion of terms within 
equation (7.2) the Reynolds equation becomes 



1 3 /o' 34>*\ 3$* f 3 /'5'„1/2 3H\ 

y.2 3Y " 2 \dX Vn 3X/ 


_3_ /£. 1/2 n 12Ub 3 (pH) 

k2 3Y Vn “ 3Y/ j Rx 3X 


(7.10) 


Relaxation methods are usually used to solve equation 
(7.10) numerically. In the relaxation process the first step is 
to replace the differentials in equation (7.10) by finite dif- 
ference approximations. The area to be considered is shown in 
Figure 7.1. The relaxation method relies on the fact that a 
function can be represented with sufficient accuracy over a 
small range by a quadratic expression. With this standard 
finite central-difference representation equation (7.10) can be 
rewritten as 


i.j 1+1. j i.j i.j-i i.j 1-1. j i.j i.j+i 


T * 

1,3 


" M? . = 0 
1,3 


(7.11) 
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where 
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Central differences are used in obtaining equation (7.11) be- 
cause they produce better estimates than forward or backward 
differences. Central differences tend to average; forward and 
backward differences tend to "overshoot" and "undershoot" the 
value of the function. More details of the derivation of equa- 
tion (7.11) are given in Hamrock (1976). 

* 

An example of a nodal structure that can be used in evalua- 
ting the pressure distribution in an el astohydrodynamic ellipti- 
cal contact is shown in Figure 7.3. A nodal structure is con- 
sidered to be suitable when the minimum film thickness within 
the conjunction aoes not change materially either when addition- 
al nodes are placed in the semimajor and semiminor axes or when 
the distances from the center of the contact to the edges of the 
computing zone are extended. Note that because of the dimen- 
sional ization of X and y the Hertzian contact area is repre- 
sented by a circle of radius unity regardless of the ellipticity 
parameter. • • 

From Figure 7.3 the features of this nodal structure are 


m 


= 4, n = 1.15, n = 1.6 


c = 5, d =■ 13 


(7.12) 


These values are used to define the nodal structure for most of 
the results presented in Chapters 8 and 9. The only exception 
is for a number of high-speed cases when the constant used to 
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detemiine the inlet distance in has to be extended from 4 to 6, 
with the other values being held constant. 

7.4 Boundary and Initial Conditions 

The following boundary conditions have been adopted: 

(1) At the edges of the rectangular zone of computation 

(Figure 7.3) the pressure is zero. This implies that <l>* is 
also zero. Specifically this means that along the bottom of 
Figure 7.3 , = o, along the left side 4>* . = 0, along the 

top and along the right side = 0* 

(2) At the cavitation boundary 

^ 3X 3Y ^ 

This condition, commonly known as the Reynolds condition, is 

'fc - 

satisfied by simply resetting <!’. v to zero whenever it occurs 

1 J J 

as a negative value. 

The pressure outside the Hertzian contact area is initially 
assumed to be zero, and therefore 0* is also equal to zero. 
That is, 3 >* = 0 when 

(X - £)2 + (Y - 1)2 1 

The pressure inside the Hertzian contact area is initially as- 
sumed to be Hertzian. Therefore from equation (3.5) the initial 
value of 4>* inside the contact area is 
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where 




* 


r 

2irabE ' 


(X - £)2 - 



(7.13) 


- 2 — 2 
(X - m) + (Y - £) <1 


7.5 Relaxation Method 

Equation (7.11) represents a system of simultaneous 
equations. There will be as many unknown values of <i)* as 
there are nodal points and equations. This system of equations 
can thus be solved by the Gauss-Seidel iterative method. In 
calculating each nodal value of this method uses the most 

recently calculated values at the surrounding nodes involved in 
the application of the finite difference equation (7.11). That 
is, nodal values calculated on that cycle are used if avail- 
able. A cycle is one complete set of calculations for all the 
nodes in the structure. The concept of the Gauss-Seidel ap- 
proach is that the best information available should be used 
instead of waiting until the next cycle. If subscript n is 
the iterant and is the particular solution to be found, 

the relaxation formula can be expressed as 
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i.j 1+1. j,n 
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i.j i.j-l.n+1 
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i, j+l,n 




* 
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(7.14) 

where is a relaxation factor. 

Therefore equation (7.14) is used, starting with node (2,2) 
and then continuing with (3,2), . . ., (IT, 2) followed by (2,3), 
(3,3), . . ., (Ti,3) and ending with (2 ,n"), (3,(7), . . ., (M,N), 
where 


M 


(m + n)d 


N = 2lc 


(7.15) 


The relaxation procedure described by equation (7.14) is contin- 
ued until 

N M 

<I>* - I 

0.1 

i.j .n+1 

j*“2, 3, . . . 3, . ♦ . 



The relaxation method provides values of o. . for every 
node within the nodal structure. Once a solution has been ob- 
tained, the dimensionless pressure can be computed from the re- 
lationship 
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With these new values of the dimensionless pressure, revised 
values of the dimensionless viscosity, density, and film thick- 
ness can be evaluated. Thus the coefficients of equation (7.14) 


(Af ., B? Cf of if and Mf J may also change. Ac- 
cordingly it is necessary to reenter the relaxation loop. This 
process is continued until the following inequality is satisfied: 


N 



j“2,3, • . • 



P - P . 

i,j,n+l x>J,n 


^i, j ,n+l 


0.1 


7.6 Normal Applied Load and Flow Rate 

The initial value of the constant Hq in equation (7.5) 
has to be estimated, but thereafter the task is to find the cor- 
rect value. To do this, the integrated normal load must be 
evaluated, where 



According to Siiripson's rule this double integral can be written 


as 



(7.18) 


where 


qi = 0 

if 

i 

is 

odd 

qi = 1 

if 

i 

is 

even 

q2 = 0 

if 

j 

is 

odd 

q2 = 1 

if 

j 

is 

even 


A new value of Hq can be estimated from the following 
expression and the current values of and F and the ap- 
plied load: 



(7.19) 


The film thickness can now be recalculated at each node for the 
new value of Hq from equation (7.5), and reentry into the 
relaxation process is required. This process is continued until 
the following convergence is satisfied: 

~ < 0.0005 

r 
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Once this convergence criterion has been satisfied, the values 
of pressure and film thickness in and around the contact are 
estoDlisnea. 

The mass flow rate per unit width in the x and y direc- 
tions shown in Figure 7.1 can be written as 


uph 


ph^ 9p 
12n 3x 


q = . £h£i£ 
y 12n 3y 

These equations can be written in dimensionless form as 


Q 


X 


UpH 


pH^ ^r\ 

12ln \b/ 9X 


(7.20) 



/ r \ d? 
12H \a) 3Y 


(7.21) 


Furthermore the total flow rates and flow angles at any node can 
be written as 




(7.22) 


Y 


a 




(7.23) 
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7.7 Flow Charts 


Figures 7.4 and 7.5 are flow charts for the numerical solu- 
tion on the digital computer of the equations developed in the 
analysis. Figure 7.4 is the flow chart for the main program. 
There are essentially three loops within the main program: the 

relaxation loop, the pressure loop, and the normal-load loop. 

In the relaxation loop , is generated. In the pres- 

sure loop the new values of . .. of the relaxation loop 

1 j J fH 1 

result in new values of pressure P. . . , , which in turn re- 

I , J jH i 

suit in new values of film thickness H. ., viscosity If. 

1 J J * » J 

and density p". The final loop, the normal-load loop, 

* » J 

ensures that the integrated normal applied load agrees with the 
initially specified value. 

Figure 7.5 is the flow chart for the film thickness subrou- 
tine SUBb. A number of calculations occur only once and need 
not be repeated on reentering this subroutine. With each new 
pressure distribution the elastic deformation is recalculated, 
and this subroutine is left with a new film thickness and there- 
fore a new value of <b* . . 

i,j,n+l 

7.8 Closure 

In this chapter a procedure for the numerical solution of 
the complete, isothermal elastohydrodynamic lubrication problem 
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for elliptical contacts has been outlined. This procedure calls 
for the simultaneous solution of the elasticity and Reynolds 
equations. In the elasticity analysis the contact zone was di- 
vided into equal rectangular areas, and it was assumed that a 
uniform pressure was applied over each area. In the numerical 
solution of the Reynolds equation the parameter a>* equals 
where P is the dimensionless pressure and H the 
dimensionless film thickness, was introduced in order to ease 
the relaxation process. The nodal structure and the boundary 
and initial conditions were given and a Gauss-Seidel relaxation 
method was used. The computer program has three major loops: 
the relaxation loop, the pressure loop, and the normal-load 
loop. The last loop requires the integrated hydrodynamic 
load-carrying capacity to be in agreement with the applied load 
within some specified tolerance. When all three loops converge, 
the pressure and film thickness in and around the elliptical 
contact area are established. 
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SYMBOLS 

constant used in equation (3.113) 
relaxation coefficients 

2 

drag area of ball, m 

semimajor axis of contact ellipse, m 

a/2m 

total conformity of bearing 
semiminor axis of contact ellipse, m 
b/2m 

dynamic load capacity, N 
drag coefficient 
constants 

19,609 N/cm^ (28,440 Ibf/in^) 

number of equal divisions of semimajor axis 

distance between race curvature centers, m 

material factor 

defined by equation (5.63) 

Deborah number 
ball diameter, m 

number of divisions in semiminor axis 
overall diameter of bearing (Figure 2.13), m 
bore diameter, m 
pitch diameter, m 

pitch diameter after dynamic effects have acted on ball, m 
inner-race diameter, m 
outer-race diameter, m 
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modulus of elasticity, Hlrtr 


E' 


1 


E 

E 

I 

e 

F 

F* 

F 

F 

F_ 


max 


f 

fb 

fr 

G 

6* 

G 

g 



1 - v! 


1 - V 


N/m‘ 


effective elastic modulus, 2 

2 2 

internal energy, m /s 
processing factor 

[("min - 

elliptic integral of second kind with modulus (1 - 1/k ) 

approximate elliptic integral of second kind 

dispersion exponent 

normal applied load, N 

normal applied load per unit length, N/m 

lubrication factor 

integrated normal applied load, N 

centrifugal force, N 

maximum normal applied load (at iJj = 0), N 

applied radial load, N 

applied thrust load, N 

normal applied load at angle N 

p 

elliptic integral of first kind with modulus (1 - 1/k ) 

approximate elliptic integral of first kind 

race conformity ratio 

rms surface finish of ball, m 

rms surface finish of race, m 

dimensionless materials parameter, aE 

2 

fluid shear modulus, N/m 
haroness factor 

2 

gravitational constant, m/s 


1/2 


1/2 
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9e 

9V 

H 

A 

H 


c,s 


min 


U 

min,r 


u 

min,s 


"min 


^min 


h 

h. 




8/3 2 

dimensionless elasticity parameter, W ' /U 

3 2 

dimensionless viscosity parameter, GW /U 
dimensionless film thickness, h/R^ 
dimensionless film thickness, H(W/U)^ = F^h/u^nQRj 
dimensionless central film thickness, h^/R^ 
dimensionless central film thickness for starved 
lubrication condition 
frictional heat, N m/s 

dimensionless minimum film thickness obtained from EHL 
elliptical-contact theory 

dimensionless minimum film thickness for a rectangular 
contact 

dimensionless minimum film thickness for starved 
lubrication condition 

dimensionless central film thickness obtained from 
least-squares fit of data 

dimensionless minimum film thickness obtained from 
least-squares fit of data 

dimensionless central-film-thickness - speed parameter, 



dimensionless minimum-film-thickness - speed parameter, 

new estimate of constant in film thickness equation 
film thickness, m 
central film thickness, m 
inlet film thickness, m 
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film thickness at point of maximum pressure, where 
dp/dx = 0, m 

minimum film thickness, m 
constant, m 

diametral interference, m 

2 

ball mass moment of inertia, m N s 
integral defined by equation (3.76) 
integral defined by equation (3.75) 
function of k defined by equation (3.8) 
mechanical equivalent of heat 

2 

polar moment of inertia, m N s 
load-deflection constant 
ellipticity parameter, a/b 
approximate ellipticity parameter 
thermal conductivity, N/s °C 
lubricant thermal conductivity, N/s °C 
fatigue life 
adjusted fatigue life 

reduced hydrodynamic lift, from equation (6.21) 
lengths defined in Figure 3.11, m 

fatigue life where 90 percent of bearing population will 
endure 

fatigue life where 50 percent of bearing population will 
endure 

bearing length, m 

constant used to determine width of side-leakage region 
moment, Nm 
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M 

M 

M 


9 

P 

s 


m 


m* 


m 

in 


N 

n 

n* 

P 

’’d 


Hz 


max 


IV, as 


m 


gyroscopic moment, Nm 

-0 75 

dimensionless load-speed parameter, WU 
torque required to produce spin, N m 
mass of ball, N s^/m 

dimensionless inlet distance at boundary between fully 
flooded and starved conditions 
dimensionless inlet distance (Figures 7.1 and 9.1) 
number of divisions of semimajor or semiminor axis 
dimensionless inlet distance boundary as obtained from 
Wedeven, et al. (1971) 
rotational speed, rpm 
number of balls 
refractive index 

constant used to determine length of outlet region 

dimensionless pressure 

dimensionless pressure difference 

diametral clearance, m 

free endplay, m 

2 

dimensionless Hertzian pressure, N/m 

2 

pressure, N/m 

2 

maximum pressure within contact, 3F/2uab, N/m 

2 

isoviscous asymptotic pressure, N/m 
solution to homogeneous Reynolds equation 
thermal loading parameter 

2 

dimensionless mass flow rate per unit width, qiiQ/pQE'R 
reduced pressure parameter 

2 

volume flow rate per unit widtn in x direction, m /s 
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R 

R. 


r r 
ax’ ^bx’ 

""ay* *^by , 

''c- ^ 

■■s- 

r 

S 

S* 

S 


0 


s 

T 

T 

Tm 

^b 

Tf* 
at* 
T, 


1 


t 

t. 


2 

volume flow rate per unit width in y direction, m /s 
curvature sum, m 

arithmetical mean deviation defined in equation (4.1), m 
operational hardness of bearing material 
effective radius in x direction, m 
effective radius in y direction, m 
race curvature radius, m 

radii of curvature, m 

cylindrical polar coordinates 

spherical polar coordinates 

defined in Figure 5.4 

geometric separation, m 

geometric separation for line contact, m 

empirical constant 

shoulder height, m 

^O^l^max 

tangential (traction) force, N 
temperature, "C 
ball surface temperature, “C 
average lubricant temperature, “C 
ball surface temperature rise, °C 

viscous drag force, N 
time, s 

auxiliary parameter 

velocity of ball-race contact, m/s 
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velocity of ball center, m/s 

dimensionless speed parameter, tiqU/E'Rj^ 

surface velocity in direction of motion, (u^ + U(j)/2, m/s 

number of stress cycles per revolution 

sliding velocity, u^^ - U|j, m/s 

surface velocity in transverse direction, m/s 

dimensionless load parameter, F/E'R^ 

surface velocity in direction of film, m/s 

dimensionless coordinate, x/R^ 

dimensionless coordinate, y/R^ 

dimensionless grouping from equation (6.14) 

external forces, N 

constant defined by equation (3.48) 

viscosity pressure index, a dimensionless constant 

coordinate system 

. . 2 

pressure-viscosity coefficient of lubrication, m /N 
radius ratio, R„/R„ 

J ^ 

contact angle, rad 

free or initial contact angle, rad 

iterated value of contact angle, rad 

curvature difference 

viscous dissipation, N/m s 

total strain rate, s"^ 

elastic strain rate, s”^ 

viscous strain rate, s“^ 
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6 

6 * 


6 


0 


6 


r 


6 


t 


6 


X 


6 


6 


e 


^2 


T1 

n 


”0 

’Ico 

e 

A 


X 


flow angle, deg 
total elastic deformation, m 
lubricant viscosity temperature coefficient, 
elastic deformation due to pressure difference, m 
radial displacement, m 
axial displacement, m 
displacement at some location x, m 
approximate elastic deformation, m 
elastic deformation of rectangular area, m 
coefficient of determination 
strain in axial direction 
strain in transverse direction 
angle between ball rotational axis and bearing 
centerline (Figure 3.10) 
probability of survival 

absolute viscosity at gauge pressure, N s/m^ 
dimensionless viscosity, ti/tiq 

O 

viscosity at atmospheric pressure, N s/m 
6.31x10"^ N s/m2(0.0631 cP) 
angle used to define shoulder height 
film parameter (ratio of film thickness to composite 
surface roughness) 

equals 1 for outer-race control ana 0 for inner-race 
control 

second coefficient of viscosity 

Archard-Cowking side-leakage factor, (1 + 2/3 a^)“^ 

relaxation factor 
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V 


5 


p 

"p 

^0 

a 


T 





’I' 

«0 

(I) 

"b 

"b 


coefficient of sliding friction 
7/n 

Poisson’s ratio 

divergence of velocity vector, (au/ax) + (av/ay) + (aw/az), s“^ 
lubricant density, N s^/m^ 

dimensionless density, p/pq 

. ? 4 

density at atmospheric pressure, N s‘-lm 

p 

normal stress, Hlmr 

2 

stress in axial direction, N/m 

p 

shear stress, N/m 

2 

maximum subsurface shear stress, N/m 
shear stress, N/m^ 

2 

equivalent stress, N/m 

2 

limiting shear stress, N/m 

ratio of depth of maximum shear stress to semiminor axis of 
contact ellipse 
Ph3/2 

<»>k=l 

auxiliary angle 
thermal reduction factor 
angular location 
limiting value of f 

absolute angular velocity of inner race, rad/s 
absolute angular velocity of outer race, rad/s 
angular velocity, rad/s 

angular velocity of ball-race contact, rad/s 
angular velocity of ball about its own center, rad/s 
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Subscripts: 

a 

b 


c 

be 

IE 

IR 

i 

K 

min 

n 

0 

PVE 

PVR 

r 

s 

x,y,z 

Superscript: 

(“) 


angular velocity of ball around shaft center, rad/s 
ball spin rotational velocity, rad/s 

solid a 
solid b 
central 
ball center 

isoviscous-elastic regime 

isovi scous-rigid regime 

inner race 

Kapitza 

minimum 

iteration 

outer race 

piezoviscous-elastic regime 
piezovi scous-rigid regime 
for rectangular area 
for starved conditions 
coordinate system 

approximate 
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Figure 7. 2. - Components of film thickness for 
elliposidal solid near plane. 







Figure 7.4. - Flow chart of main program. 












Figure 7 . 5 . * Flow chart of subroutine SUB6. 
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